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Abstract: 

We propose a systematic approximation scheme for solving GLAP evolution equations at 
small Bjorken-x. The approximate solutions interpolate smoothly between hard (singular 
as aW 1+A ), A > 0) and soft (singular at most as x _1 ) initial conditions and may be applied 
in a wide range of Q 2 . The small-x behavior of F$(x, Q 2 ) which is extracted from a fit to 
HERA data agrees with results from a global fit. 
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Understanding the small-x HERA data [1|, |2| is nowadays one of the most impor¬ 
tant challenges to perturbative QCD In a conventional approach the link between 
theory and the data is provided by a mathematical parameterization constructed to obey 
QCD evolution according to Gribov-Lipatov-Altarelli-Parisi (GLAP) equations. Although 
global fits in the whole x region play necessarily the most important role B I- i 0 , many 
approximations aimed specifically at the small-x region I, B 0 0 have been proposed 
with a goal to provide more insight |H| into the otherwise numerically complicated analy¬ 
sis. As it has been emphasized in [|13j such an approximate procedure can be justified only 
if its accuracy with respect to the exact calculation is higher than a typical accuracy of 
the experimental data. In this letter we present an analysis based on a standard analytical 
approach to the LO GLAP evolution equations which includes corrections necessary for 
satisfactory agreement with the state-of-art numerical integration for both hard, (singu¬ 
lar as x - ( 1+A ), A > 0) and soft (singular at most as x -1 ) initial conditions. By fitting 
the final formula to the data we can explicitly demonstrate that the analytical approach 
predicts the same small-x shape of the input distributions as the global fit, as it should 
be. Our analysis also shows that at Q 2 evolution lengths which are typical for the present 
data, hard initial conditions may produce significant contributions which possess a simi¬ 
lar small-x behavior as those resulting from soft initial conditions. As a consequence, we 
have found that the same data can be described by either moderately hard A ~ 0.2 or 
soft initial conditions, depending on the starting point of the evolution. 

Let us now describe the analytic approach to the LO GLAP evolution which we apply 
later to the analysis of the small-x HERA data. As usual, we denote as 


Qs(n,t ) = 

g(n,t) = 


dxx n q s (x , t ), 
dx x n g(x, t ) 


( 1 ) 


the moments of the quark-singlet q s (x,t ) and gluon g(x,t) distributions at a scale t = 
ln(Q 2 /Aq CD ). At leading order the GLAP evolution of the quark singlet distribution in 
terms of the initial distributions q s {n, to) and g(n,to) has the well-known form WA 


q s (n , t) = {(1 - h 2 (n)) q s (n, t 0 ) - hi(n) g(n , t 0 )} e^o A+(n)c 

+ {h 2 (n) q s {n , t 0 ) + hi(n) g(n, f 0 )} e^ A ~ (n)C . 


( 2 ) 


In the above expression ( = ln(a(Ql)/a(Q 2 )) denotes the Q 2 evolution length and A ± 
are the eigenvalues of the anomalous dimension matrix. The coefficients hi and h 2 may 
be expressed through A± and the anomalous dimensions as 


hi(n) = 


A'qgi 


n 


h 2 (n) = 


7 qqH ~ A +(n) 


( 3 ) 


To obtain the Bjorken-x distribution from 
which reads 

. , 1 

xq s {x,t) = 


A_(n) — A + (n) ’ A_(n) — A + (u)’ 

one has to take the inverse Mellin transform 


rc-\-ioo 


2ni Jc-i 


dnx n q s (n,t), 


( 4 ) 
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where the integration contour runs to the right of all singularities of q s (n,t). As we shall 
be ultimately interested in the small-x behavior of the quark distribution, let us assume 
that the initial conditions are given by simple poles 


qs(n,t 0 ) 


Mto) 

n - X q ’ 


9{n, to) 


Ag(t 0 ) 

n-Xg 1 


( 5 ) 


with 1 > X q g > 0. Anticipating that such a power-like rise in the initial condition is only 
valid for small enough x, one could introduce an additional parameter Xo above which 
one neglects the parton distribution and below which one assumes a power-like behavior. 
With this interpretation in mind it is evident that 0 < Xo < 1. This yields the following 
modified a n sat z for the initial boundary, 


Qs(n, t 0 ) 


Mto)x 0 ,g q 
n - X, “•« ’ 


g(n,t 0 ) 


Ag(to)xo X g g 

n-Xg 


( 6 ) 


where we have introduced separate x 0 ’s for quark and gluon distributions. As far as the 
magnitude of x 0 ’s is concerned Ref. || convincingly argues that x 0 , g should be of order 
of 0.1. 

We have found that the introduction of x 0g , g causes problems. From a practical point of 
view the model parameters become highly correlated (cf. equation (^)). Secondly, the 
HERA data S! favor Xo, q to be larger than one which contradicts the interpretation of 
xo as the boundary of the power-like behaviour. For these reasons we have decided to use 
instead the following ansatz for the initial conditions 


q 8 (n,t Q ) = Aq ^ + B q (t 0 ) , g(n,t Q ) = + B g (t 0 ). (7) 

The constants B q g model the leading corrections to a power-like small-x boundary. Fur¬ 
ther corrections nB^(to) + ■ ■ ■ can be introduced if necessary. 

Our strategy for solving the Mcllin inversion (|4|) with initial conditions (|7]) will be the 
usual one, i.e., first we close the contour in the left-half plane and then successively the 
residues are taken, starting from those with the largest real parts. The smaller x is, 
the less residues will be needed. Therefore we first discuss the poles inside the unit disc 
n < 1. Apart from the pole at n = A the only other pole singularity inside the unit disc 
occurs at n — 0. This can be read off from the Laurent expansions of A± and h\^ which 
start like 


ld+ (n) 

A» 
A_(n) 
hi(n), 1 - h 2 (n) 


- + A », 

n 

- 6+ + 0(n ), 

^(l), 

O(n), 


( 8 ) 
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where a = 12 //3 0 and <5+ = (11 + 2nj /2T)/ /3 0 . The functions A^, A_ and h]/ 2 are analytic 
inside the unit disc. From this one concludes that the contribution to the Mcllin integral 
([]]) from the A_-term in (| 2 ;) is just the residue of the simple pole at n = X 

xq.(x, Q 2 ) = A q (t 0 ) h 2 (X q ) e a ' A -( x i)+ bX Q + Ag (t 0 ) hi(X g ) e a ' A -( A 9)+ bA s ; ( 9 ) 

where a' = 2(//3q and b = In (1). The Mellin inversion of the A + -term in (D) is complicated 
through the additional essential singularity at n = 0. Using that A^, and h\ i2 are analytic 
inside the unit disc one can formally rewrite the Mcllin inversion of the A + -term as 

xq + (x,Q 2 ) = A q (t 0 )U q (d b ,() J(a,b;\ q ) + A g (t 0 )Ug(d b ,() J{a,b;X g ), 

+ {B q (to)U q (d b ,() + B g (t 0 )Ug(d b ,()} (^) 2 h(V^b), 

1 r en +bn 

J{a,b]X) = — / dn -- (10) 

2m J'y 7i — X 

where again a = 12C/An b = ln(-), I\ is the modified Bessel function |15|, and 7 
surrounds both poles with positive orientation. The differential operators in (JT^) are 
constructed in such a way that their action reproduces the expansion (||), e.g., 

7(a,0 = -Ma7 A »>, ( 11 ) 

where d b — ^ denotes the derivative over the parameter b. Note, that we do not expand 
the initial condition around n = 0. Such an expansion would result in a term proportional 
to 1/A which diverges at A = 0. It is clear that to obtain a sensible soft initial condition 
limit one has to keep the initial condition term as a whole. 

Now, from properties of modified Bessel functions it can be shown that J(a, 6 ; A) 
admits the following representation 



In the limit ( = 0 = a, i.e. no Q 2 evolution, the above formula exactly reproduces the 
corresponding initial condition. 

As the next step consider the operator U(d b , (). For fixed ( the function U(y, () should 
be looked at as a Taylor series in y. After substituting y by d bl each term y l in the above 
Taylor expansion becomes a differential operator that acts on J(a, 6; A) to give 

1 r . e™ +bn 

Ji(a,b- A) = — / dnn 1 --. (13) 

I'm j'y 71 — X 

Again, the above integral can be evaluated to give 


Ma,b;X) = X l J2 

k =0 


t) 


Ik{2y/ab) + 


A i+1 
2 a 


,b\ 


r2x/~ab 


A 2 

duue~^ u I 0 (u) 


(14) 
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( 15 ) 


where n 0 = </|. For a soft boundary Jj(a, 6 ; A = 0) simplifies to 


Jj(a, 6 ; A = 0) = n l 0 Ii(2\/ab). 


Combining the representation 


h(z) = ~ I dt e zcost cos (kt) 

71 JO 


(16) 


with the Taylor expansion of U(db,() if is possible first to compute the geometric sums 
over i in (0) and then to evaluate their infinite series by a single complex integral. One 
obtains 


xq+(x,t) = A q (t 0 ) \ q U q (\ q ,()-^T L duue I 0 (u) 


+ Ag(t Q ) \gUg(\g,C) 


gb\q i‘2\/ab 

2 a Jo 

£>b\ g ~2\Jab 


3 . „2 


2 a Jo 


duue 4 “ “ Io(u) 


— [* dt e 2 '^ >cost Re ( ^g C) A g(jo) UgiJ^gi C) 

n Jo l Ag — n Q e u A„ — n 0 e i1: 


^0 r dt e 2Vabcost Re f A M Uqjnoe*, 0 e u + A g (t 0 ) U g {riQe it X)e i 


7T JO 

no 


A q - 7l 0 e i 


A q - n 0 e i 


it 


+ “ Jq dt e 2 ' A ^ cost Re ^B g (t 0 ) U q (n 0 e it , () + B g (t 0 ) U g (n 0 e u , C)} , (17) 

where Re denotes the real part. This formula is valid if U(z, Q admits a Taylor expansion 
around z = 0 along the integration contours z = n o exp (it). Thus, because the conver¬ 
gence radius of Ug !q (z, () is the unit disc, only data with no < 1 should be compared with 
( pT|) which is not a severe restriction in practice . Of course the same restriction also 


applies to approximations that are constructed from (h4|) , i.e. where only a finite number 
of terms are taken into account. The formula ([17]) looks more complicated than it really 
is, although it is necessary to take care about a proper definition of the multi-valued 
complex square root used in the definition of the anomalous dimension along the contour 
of integration. In the soft boundary limit the first two terms in ([17]) vanish and the last 
term simplifies somewhat. Apart from the kinematical restriction tt-o < 1, equation 0 ) 
together with (j!|) yield the exact contribution to the parton radiation arising from the 
leading singularity at n = 0 of the anomalous dimensions and the leading pole at n = A 
of the boundary. 

The next singularity in the Mellin plane at n — — 1 can be treated along the same lines 
that lead to dnp but gives only a numerically small correction at small x because of an 
extra factor x that multiplies most terms. An important difference now is the absence of 
a simple pole like the one at n = A above. It can explicitly be shown that this causes all 
contributions from the n — — 1 singularity to vanish in the no-evolution limit as it should, 
because the initial boundary is already recovered by the leading poles at n = A and n — 0 . 
Therefore one can expect that for moderate Q 2 this term is additionally suppressed. In 
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principle it is an easy task to incorporate effects from the pole at n = — 1 but for con¬ 
sistency one should then also include the next pole from the initial conditions. Rather 
than to proceed this way we neglect all poles with Re{n ) < —1 and check the quality of 
the approximation a posteriori by comparing with the exact result. Figure 1 shows such 
a comparison with the exact LO CTEQ |J at Q 2 = 10 GeV. The approximation (jTj) to 
the initial conditions has been determined from the CTEQ input parameterizations. For 
g_ we use equation (|9j) and q + is constructed by taking in ([T4]) the first non-vanishing 
term in the expansion of U q:g (y,() around y = 0. One learns from this figure that when 
keeping only data with, say, x < 5 • 1CT 3 , it is certainly justified to neglect poles with 
Re(n) < —1. Further, despite the fact that typically no is not particularly small [|13j, it is 
not necessary to evaluate the resummed formula (JT7|) but it is sufficient to take just the 
first term from the expansions that are generated by (0). u. 


xq s (x,t) = 


(0.198 A q (t 0 ) + 0.444 A 5 (f 0 )) (Jj 2 / 1 (2v / ^) 
(0.198 X q A q (t 0 ) + 0.444 A 9 A g (t 0 )) I 0 (2y/ab) 


(0.198 B q (to) + 0.444 B g (t Q )) J l2{2\fab) 

A 2 , fZVab Xq 2 
0.198 A q (to) x q duue 4a u Io(u) 
2 a Jo 

A? , r^Vab \ 2 

0.444 A g (t 0 ) —^ x~ 9 / duue~^ u Iq(u) 

2 a Jo 

’(1 - 0.198AQ A q (t 0 ) x~ Xq - 0.444 X g A g (t 0 ) x 


= -5+C 


,- 6 +1 


r. - 


^l.l 85C 


(18) 


This is largely due to the fact that the Taylor coefficients that multiply these terms are 
rather small for sub-leading terms and values of ( typical for the current data. Note, 
that the extreme case where only the first non-vanishing term in the expansion of U qt9 is 
kept and X g , X q are set to zero corresponds to the soft boundary solution that has been 
obtained in || by means of a saddle-point approximation for large ( and large £. 

Furthermore we compare in Figure 1 our expansion scheme with a naive pole approx¬ 
imation which takes into account the residue of the rightmost pole in Mellin space at 
n = X q , X g only. A priori we expect such an expansion scheme to work in the regime of 
large A, where the pole is sufficiently separated from the essential singularity at n = 0 aris¬ 
ing from the anomalous dimensions. Already for A values around 0.2 — 0.3 as represented 
by the CTEQ parameters (see Table [L|) the difference between exact and approximated 
solution becomes larger than the experimental errors. In case of yet smaller A we have 
checked that the discrepancy strongly increases. Thus the naive pole approximation has 


to be replaced by the more accurate formula, like equation (|lq) , in a data analysis where 
A is a variable fit parameter. 

Using the above considerations we have performed a fit |T 6 | to the HERA data [|I|. 
|2] with the formula (|T 8 |), starting from scales Qq = 1 and 2.56 GeV 2 , the latter scale 


corresponding to the CTEQ starting point. In both fits the number of quark flavors was 


5 













globally assumed to be nj = 4. The results are shown in Table [I]. For orientation, the 
table contains also the corresponding fit parameters determined from the CTEQ initial 
conditions. As far as the fit with the starting point Qq = 2.56 GeV 2 is considered, the 
parameters A q and A g as well as the slopes X g and X q agree well with what has been found 
in from a global fit. It demonstrates the self-consistency of the GLAP evolution in the 
small-x region: the parameters of the fit which determine the leading small-x behavior 
are indeed determined mainly by the small -x data. At the same time if one is interested 
only in the QCD description of the small-a; data it is legitimate to use the analytical 
approximation to the GLAP evolution, saving much effort as compared with a global 
fitting procedure. The parameters B q and B g are reproduced less accurately, but still the 
present fit seems to reproduce their correct signs. Note, that B q comes out from both the 
CTEQ and our fits greater than zero which contradicts the interpretation xo, q , as we have 
discussed above. The relatively large errors are due to the fact that B q (to) and B g (to) 
generate non-leading contributions at small-x and therefore are less sensitive to the data 
in this region. 

The fit with a low starting point of Qg = 1 GeV 2 results with the same data being 
parameterized by a soft gluon boundary, although the corresponding X g parameter is 
hardly constrained by the fit. It means that for a sufficiently large evolution length 
soft and hard boundaries result in a similar shape of F 2 (x, Q 2 ) at small-x. Note, that 
although equation (|l^) describes asymptotic behavior of the evolution starting from a 
hard boundary, the first three terms on the right-hand side have a non power-like, i.e. 
soft, small-£ behavior. To illustrate this point we show in Figure 2 various contributions 
to F 2 (x,Q 2 ) for Q 2 = 10 and 1000 GeV 2 , respectively, as determined by our fit with 


the starting scale Ql = 2.56 GeV 2 . As already stated in |K| it is important to include 


q-(x,Q 2 ) in an analysis including data in a moderate Q 2 range. This contribution is 


represented by the last term in Eq. ([l 8 |) . For large Q 2 this part of the solution gets 
strongly suppressed as expected from the very beginning. On the other hand the soft 
part, identified as the first three terms on the right-hand side of equation ([18]), has a 
shape which is similar to the full solution. In other words, in the x domain covered by the 
current experiments soft and power-like initial conditions can be indistinguishable. This 
observation may help to understand the known fact that the global fits can explain the 
data equally well starting from soft [ 6 fl or hard [|7j] initial conditions. 

Keeping X g and X q equal to zero, and allowing for the starting point of the evolution 
to be as low as 0.45 GeV 2 , we have obtained the fit presented in the last row in Table |T]. 
It is consistent with results of || where the LO fit starting from a yet lower scale 0.24 
GeV 2 and valence-like distributions with X g , X q < 0 was presented, but one should keep 
in mind that the approximation ( |T 8 | ) does not describe the GLAP evolution from such 
initial conditions with a sufficient accuracy. 

Finally we have investigated the double asymptotic scaling behavior of F 2 as function 
of the scaling variables p = \J £/£, a = v^C (£ — ln(l/x), ( = In “ s ^q °] ), as advocated 
p. After rescaling F 2 by the factor 


m 


Rf = Nf p v^exp (—2 7 cr + 8 + a/p) 


(19) 


6 





Data range 

Qo 

Ag(t 0 ) / A q (t 0 ) 

Bg(t 0 )/Bg(t 0 ) 

^g/^q 

X 2 / d.o.f. 

x < 0.005, 

Q 2 > 3GeV 2 

2.56 GeV 2 
(fixed) 

0.84T0.08, 

0.57T0.03 

-0.81il.02, 
1.21 ±2.33 

0.27±0.02, 

0.18±0.01 

0.78 

CTEQ § 

2.56 GeV 2 
(fixed) 

0.854, 

0.502 

-2.693, 

1.799 

0.305, 

0.200 


x < 0.005, 

Q 2 > 3GeV 2 

1.0 GeV 2 
(fixed) 

0.19T0.05, 

0.70T0.03 

-1.62±0.49, 
6.27 ±0.77 

0.32±0.03, 

0 .11±0.01 

0.78 

x < 0.005, 

Q 2 > 3GeV 2 

0.45 ± 0.01 GeV 2 

0.00T0.06, 

1.04T0.04 

0.00±0.17, 
1.61 ±0.29 

0.00 (fixed), 
0.00 (fixed) 

0.87 


Table 1: Results of the fit of equation (|l^) to the small-a; HERA data Si- For orienta¬ 
tion, the table contains also the corresponding fit parameters determined from the results 
of the CTEQ LO global fit |J| . In the last row we have kept X g and X q equal to zero and 
allowed for perhaps unphysically low starting point Qq in order to present a fit with soft 
initial conditions for both quark and gluon distributions. All errors are parabolic, see |]TB 
for details. 


where h + = (11 + Arif)//3 0 , we observe in Figure 3 that the data accumulate in a flat 
area in the p-cr-plane, despite the fact that the plot corresponds to a power-like initial 
quark distribution at Qq = 1 GeV 2 . In other words, the kinematic bounds of current 
experiments restrict the position of the data points to a domain in the p-cr-plane, in 
which RfF 2 hardly differs from a flat plane. On the other hand, the double-asymptotic 
scaling interpretation does not hold anymore if the starting point is taken to be Q q = 2.56 
GeV 2 |H]], mainly due to the redefinition of the Q 2 evolution length £, see Figure 4. 


In summary, we have presented an approximation scheme for solving analytically the 
GLAP evolution at small x that treats corrections to the strict small-re limit in a systematic 
way. In contrast to saddle-point approximations there is no need to have Q 2 very large. 
Our solution interpolates smoothly between soft and hard boundaries. For moderately 
hard boundaries we have found that the pole approximation alone is not sufficient to 
adequately reproduce the full QCD evolution. The fit shows that it is not necessary 
to perform a global analysis in the full Bjorken-x range if one is interested only in the 
parameters which determine the small -x behavior of parton distributions at the initial 
scale. We have confirmed the wisdom arising from the global fits, that the slopes of initial 
quark and gluon distributions depend critically on the starting point of the evolution. 
The same HERA data can be described by either a hard boundary at Q 2 = 2.56 GeV 2 or 
a soft one if the starting point is lowered below Q 2 = 1.0 GeV 2 . The interpretation of the 
current data in terms of double asymptotic scaling crucially depends on the magnitude 
of the starting point of the Q 2 evolution. 
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Fig. 1 Comparison of the full (solid line) and the approximate evolution according to 
equation (|T8|) (dot-dashed line) at Q 2 = 10 GeV 2 . The initial conditions at Q 2 = 
2.56 GeV 2 correspond to the CTEQ4L fit ||. The dashed line shows the result 
of the naive pole approximation, based on taking into account the residue of the 
rightmost pole in the Mellin plane only. Data points Jl]| correspond to Q 2 = 10 
GeV 2 . 
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Fig. 2 Contributions of soft terms, and small eigenvalue 5/18 q~(x, Q 2 ) terms to F 2 (x, Q 2 ) 
at Q 2 = 10 GeV 2 (upper plot) and Q 2 = 1000 GeV 2 (lower plot), as obtained from 
our fit with the starting scale Ql = 2.56 GeV 2 . Full lines represent F 2 (x,Q 2 ). Dot- 
dashed lines show the soft contribution, as given by the first three terms on the right- 
hand side of equation (p~8|) . Dashed lines show the small eigenvalue contribution 
5/18 q-(x,Q 2 ). 
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Fig. 3 RpF 2 as a function of the scaling variables p and a as given by the fit with 
Ql = 1.0 GeV 2 . Dots denote data points JT], |2| with Q 2 > 3 GeV 2 . Note that the 
small-x region corresponds to large values of p. 
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Fig. 4 RpF 2 as a function of the scaling variables p and a as given by the fit with 
Ql = 2.56 GeV 2 . Dots again denote data points |1], jj2j] with Q 2 > 3 GeV 2 . The 
RfF 2 (p , cr) surface is not flat and the double-asymptotic scaling interpretation does 
not hold anymore. 
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